We are IntechOpen, 
the world’s leading publisher of 


Open Access books 
Built by scientists, for scientists 


5.300 130,000 150M 


ailable International authors and editors Downloads 


Our author among the 


154 TOP 1% 12.2% 


Countries delivered to most cited s Contributors from top 500 universities 


Selection of our books indexed in the Book Citation Index 
in Web of Science™ Core Collection (BKCI) 


Interested in publishing with us? 
Contact book.department@intechopen.com 


Numbers displayed above are based on latest data collected. 
F information visit www.intechopen.com 


ay 


15 


Fourier Factorization in the Plane Wave 
Expansion Method in Modeling Photonic Crystals 


Roman Antos! and Martin Veist? 

‘Institute of Physics, Faculty of Mathematics and Physics, Charles University in Prague 
2 Institute of Biophysics and Informatics, First Faculty of Medicine, Charles University in 
Prague 

Czech Republic 


1. Introduction 


Photonic crystals are modern artificially designed periodical systems capable to affect 
the motion of photons in a similar way that the periodicity of the atomic potential 
in a semiconductor crystal affects the motion of electrons. The physical properties of 
light in a photonic crystal resemble those of electrons in atomic crystals, leading to 
forbidden propagation of electromagnetic modes at certain frequencies, as demonstrated by 
Yablonovitch (1987), Ho et al. (1990), Joannopoulos et al. (1995; 1997), etc. The existence of the 
optical band gap (which is the part of the spectrum for which the wave propagation is not 
possible) makes photonic crystals broadly interesting from many viewpoints of fundamental 
research and applications. Recent studies are motivated by promising applications such as 
purely optical integrated circuits [e.g., by Lin et al. (1998), Noda (2006), or Hugonin et al. 
(2007)], artificial metamaterials with high tunability [Datta et al. (1993); Genereux et al. (2001); 
Krokhin et al. (2002); Reyes et al. (2005)], high-sensitivity photonic biosensors [Block et al. 
(2008); Skivesen et al. (2007)], or devices based on phenomena not accessible in conventional 
media [Benisty (2009); Kosaka et al. (1998); Krokhin & Reyes (2004)]. It is also necessary for 
accounting for structural colors of wings of butterflies or beetles, feathers of birds, or iridescent 
plants [Kinoshita & Yoshioka (2005); Vukusic & Sambles (2003)]. 


Since the optical properties of photonic crystals strongly depend on their geometrical 
structure, used materials, etc., their proper design is crucial for the correct device functionality. 
Thorough theoretical analysis therefore takes place in the development, using various 
numerical methods for calculation including methods of finite difference in the time or 
frequency domain or finite element methods [Joannopoulos et al. (1995)]. One of the most 
common calculation techniques applied to photonic crystals is the plane wave expansion 
method, which is a frequency-domain approach based on the expansion of the fields and 
material parameters into the Fourier (or reciprocal) space. The components of this expansion 
represent the definite-frequency states. After some necessary truncation of the complete 
basis (plane waves with a finite cutoff), the partial differential equations are then solved as 
a linear-algebraic problem. However, the convergence rate of this method strongly depends 
on the implementation of Maxwell’s equations in the truncated plane-wave basis [Meade 
et al. (1993); Sozuer et al. (1992)]. In the case of periodic discontinuities (typical for photonic 
crystals) the convergence is rather poor so that the computer calculations might become 
extremely time- and memory-consuming. 
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Because the underlying physical phenomenon in the optical behavior of photonic crystals is 
based on diffraction (therefore the lattice constant of a periodic structure has to be in the same 
length-scale as half the wavelength of the electromagnetic wave), several conclusions can be 
advantageously adopted from the classical coupled-wave theory, one of the most effective 
methods for modeling diffraction of electromagnetic waves by periodic gratings, which was 
developed during several past decades. This can provide a high enhancement to the plane 
wave expansion method, resulting in the reduction of the computation resources. 


In the mid 1990s Li (1996) showed that Laurent’s “direct” rule, which had always 
been adopted in conventional formulations to factorize the truncated Fourier series that 
corresponds to products of two periodic functions, presents bad convergence when the two 
functions of the product are simultaneously discontinuous. He suggested three Fourier 
factorization rules (briefly summarized in Section 2.3) and applied them to one-dimensional 
(1D) diffraction gratings. This major breakthrough in the grating theory (called “fast Fourier 
factorization”) was soon applied by many authors to various grating structures with arbitrary 
periodic reliefs, anisotropic [Li (1998)] and slanted [Chernov et al. (2001)] periodic systems, 
their various combinations [Li (2003); Watanabe (2002); Watanabe et al. (2002)], and other 
systems [Bonod et al. (2005a;b); Boyer et al. (2004)]. 


Later Li (1997) applied the factorization rules to two-dimensional (2D) periodic structures 
treated by “zigzag” Fourier expansion, which yielded an improvement for rectangular dots 
or holes. However, Popov & Neviere (2000) have pointed out that the staircase approximation 
(of the coupled wave theory using the slicing of relief profiles) in combination with the 
traditional formulation of differential equation within one slice violates Li’s factorization 
rules. This was a major complication for the analysis of the periodic systems made of rounded 
elements. Therefore, they applied a coordinate transform to treat individually the normal and 
tangential components of the electric field on 1D sinusoidal-relief gratings, which enabled the 
application of the correct rule for each field component and thus improved the convergence. 


Later David et al. (2006) utilized the normal-tangential field separation to 2D photonic crystals 
composed of circular or elliptical holes. Similarly, Schuster et al. (2007) applied this method to 
2D gratings, and also suggested more general distributions of polarization bases [Gotz et al. 
(2008)]. These approaches, always dealing with linear polarizations, enabled a significant 
improvement of the convergence properties, but ignored the fact that the transformation 
matrix between the Cartesian and the normal—tangential component bases of polarization 
became discontinuous at the center and along the boundaries of the periodic cell, which 
slowed down the resulting convergence. To overcome these discontinuities, a distribution 
of more complex (i.e., generally elliptic) polarization bases was recently suggested to improve 
optical simulations of 2D gratings and photonic crystals [Antos (2009); Antos & Veis (2010)]. 


Our chapter will describe in detail the application of the Fourier factorization rules to the plane 
wave expansion method for numerical analysis of general photonic crystals. Section 2 will 
introduce the principle of the plane wave expansion together with the notation of matrices 
and factorization theorems. Section 3 will refer to 1D photonic structures made as periodic 
stratified media. The consistency of the correct factorization rules with classical theory of Yeh 
et al. (1977) and Yariv & Yeh (1977) will be shown, pointing to the correct boundary conditions 
of the tangential components of the electric and magnetic field on multilayer interfaces. 
Section 4 will repeat our previously described methodology for 2D photonic crystals made 
of circular elements, and Section 5 will generalize it to elements of other shapes. Sections 6 
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and 7 will propose how to factorize anisotropic and three-dimensional (3D) photonic crystals, 
respectively. 


2. Plane wave expansion 


2.1 General remarks 


The modes of photonic crystals are in principle the eigensolutions of the wave equation with 
an inhomogeneous, periodic relative permittivity e(r). One possible version of the wave 
equation is the equation for the unknown electric field with an unknown frequency, 


2 
Iy x (V x 2) = 5E, (1) 
where w/c? is its a (w is the frequency and c the light velocity in vacuum) and 
E(r —e ik-r D Cn ye tol (mpx-+ngy+lsz) _ =y Cn —iko(pmx+qny+s)Z) (2) 
m,n,l m,n,l 


is its eigenfunction, which has the form of a pseudoperiodic Floquet-Bloch function. Here 
p = 27/Ax, q = 27/Ay, and s = 27/A; are the normalized reciprocal lattice vectors. For 
simplicity we assume the periods A; along the Cartesian axes throughout this chapter. For 
brevity we have also defined k = ko[po, qo, S0], Pm = po + mp, Qn = qo + ng, and s; = so + Is. 
(Analogously we could write the wave equation for an unknown magnetic field H or any 
other field from Maxwell’s equations.) 


Owing to the periodicity of the problem, the plane wave expansion method is the reference 
method for the mode calculation. It is based on the Fourier expansion of the field such as in 
Equation 2 and on the Fourier expansion of a material function, either the permittivity or the 
impermittivity y(r) = 1/e(r), 


e(r) = iJ T a (3) 
m,n,l 
E Nmn eto (mpx+nqy+lsz) (4) 


The rules for choosing the most spe material parameter and the most appropriate 
field for the Fourier expansion are governed by various methods of Fourier factorization. In 
the past, the E method (7 and the electric displacement D were expanded), H method (7 and 
H were expanded), and Ho method (e and E were expanded) were the typical choices. 


2.2 Matrix notation 


Now we carry out the transformation of the partial differential equations into matrix equations 
in order to solve the eigenproblem by linear-algebraic methods. For simplicity we limit 
ourselves to 1D and 2D photonic crystals, and always choose the direction of propagation 
in the xy plane, so that 0, = 0. With these restrictions we write 


e(x, y) = 5 Emn e` i(mpx-tngy) (5) 
m,n=—0o 

feyn= >. fact. (6) 
m,n=—0o 
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where f is a component of the electric field. We will now derive the matrix expressions of the 
fundamental relations 


h(x, y) = e(x,y) f(x,y), (7) 
8x(x, y) = Oxf (x,y), (8) 
8y(x%,y) = dyf (x,y), (9) 


i.e., the relations of multiplication by a function and applying partial derivatives. Assuming 
the expansions of the new functions h = Ym, n Amn eil Pmx+qny), g, = Emn &x,mn e™ilPmx+qny) 
and gy = Linn Qy,mn ei(Pm*+4n¥) we rewrite Equations 7-9 using the convolution rule, and 
applying the partial derivatives as follows: 


+00 
hmn = 3 Em-k,n-l fki» (10) 
k,l=—co 
&x,mn = —ipmfmn, (11) 
Symn = —ignfmn- (12) 


Assuming furthermore à finite number of the retained Fourier coefficients, i.e., using the 
summation Eo M T __y We can renumber all the indices to replace the couple of two 
sets m € {-M, -M +1, ..., M} andn € {-N, -N +1, ..., N} by a single set of indices 
& € {1, 2, ..., Amax}, wh &max = (2M + 1)(2N +1), alad. 


a(m,n) =m+M+1+(n+N)(2M+1), (13) 
n(a) = (a — 1)div(2M +1) - N, (14) 
m(a) = (a —1)mod(2M +1) — M, (15) 


where “div” denotes the operation of integer division and “mod” the remainder (the modulo 
operation). Then we can rewrite Equations 10-12 into the matrix relations 


[h] = lell) (16) 
[gx] = —ip [f], (17) 
[gy] = -iq [f], (18) 


where [f], [h], [gx], and [gy] are column vectors whose ath elements are the Fourier [m,n] 
elements of the functions f, h, gx, and gy, indexed by a(m,n) defined in Equation 13, and 
where |e], p, and q are matrices whose elements are defined 


Leap = &m(a)—m(B),n(a)—n(B)* (19) 
Pup = Pin(a)%eBr (20) 
Yap = In(a) ap) (21) 


where the indices on the right hand parts are defined by Equations 14 and 15 and where dy 
denotes the Kronecker delta. As a summary we can say that the multiplication by a function 
is in the reciprocal space represented by the matrix [|e] (in the sense of the limit amax — ©°) 
and that the partial derivatives are represented by the diagonal matrices —ip and —iq. 
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For 1D periodicity we choose the inhomogeneity along the y axis. In this case the a index is 
not necessary because the n index is sufficient, 


e(y) = L ene E, (22) 
fa) = »2 fae oe (23) 
lel nk = En—kr (24) 

lnk = GnOnk- (25) 


Here |e] is a Toeplitz matrix (a matrix with constant diagonals). 


2.3 Simplified theorems of Fourier factorization 


Although the theorems were derived by Li (1996) for 1D periodic functions, we here 
summarize them in the matrix formalism independent of the number of dimensions. Let f, h, 
and e be piecewise-continuous functions with the same periodicity related 


h=ef, (26) 
and let [f], [h], and [e] denote their matrices as defined in Section 2.2. 


Theorem 1. If € and f have no concurrent discontinuities, then the Laurent rule applied to 
Equation 26 converges uniformly on the whole period and hence 


[h] = fel (27) 
can be applied with fast convergence. 


Theorem 2. If £ and f have one or more concurrent discontinuities but h is continuous, then 
Equation 26 can be transformed into the case of Theorem 1, 


fain, (28) 


and hence 
Lf] = Bi [nh]. (29) 


Accordingly, we can state 


= |=] ef (30) 


referred to as the inverse rule. We say that the functions € and f have complementary 
discontinuities. 


Theorem 3. If none of the requirements of the first two theorems are satisfied, then none of the 
rules can be applied correctly because Equations 27 and 30 are no longer valid at the points 
of discontinuities, which considerably slows down the convergence. Therefore, we should 
carefully analyze the continuity of the functions and transform all the partial differential 
formulae to the first two cases. 
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3. One-dimensional photonic crystals 
3.1 Geometry of the problem 


Fig. 1 shows the geometrical configuration of a 1D photonic crystal made as periodic 
alternation of two different layers with relative permittivities €y and gy, thicknesses dı and 
dy, with periodicity A along the y axis. The relative thickness of the first layer (with respect 
to the period) is denoted w; the relative thickness of the second layer is then 1 — w. The 
coordinate system is chosen to get the uniform problem along the z axis. This means that the 
plane of incidence (plane defined by the vector of periodicity and the wave vector of incidence 
k = ko|po, 40,9], here only hypothetical since the photonic crystal is infinite) coincides with 
the xy plane. 


plane of incidence 


-iko (pPox+qny) 


-iko (Pox+qry) a,e 
1 


dj=wA ‘d,=(Il-wA % d, 
Fig. 1. Geometry of a 1D photonic crystal made as periodic alternation of two layers 


Then we distinguish between two polarizations of electromagnetic fields which can be treated 
independently. The transverse electric (TE) polarization has E perpendicular to the plane of 
incidence (Ez, Hy, and Hy are nonzero). The transverse magnetic (TM) polarization has H 
perpendicualar to the plane of incidence (Hz, Ey, and Ey are nonzero). 


3.2 Application of Fourier factorization 


Propagation in a 1D periodic medium, whose inhomogeneity along the y-axis is described 
by the relative permittivity function e(y), is governed by Maxwell’s equations (choosing the 
coordinate system uniform along the z-axis, i.e., dz = 0) 
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where the first and the second row corresponds to the TE and the TM polarization, 
respectively. The formal labels of the relative permittivity, (£) and (Z), means that the 
multiplication with the corresponding component of the electric field will be (according to 
the factorization rules) treated by either the Laurent rule (£) or the inverse rule (Z). This is 
due to the fact that Ez and Ey are continuous functions (because they are tangential to the 
discontinuities of e), whereas Ey is the component normal to the discontinuities and hence 
discontinuous (while the product £E; is continuous). By separating one field in each case, we 
obtain the wave equations 


= (a + a2) Ez = Nees, (TE) (33) 
1 = z 
z (Le + 8y1ay ) H, =H, (TM) (34) 


or, after expanding the permittivity or impermittivity and the fields into the Fourier series, the 
matrix formulae 


lel? (p2 +q?) [Ez] = S [Ez], (TE) (35) 
([2] 73+ alela) [He] = SH]. (TM) (36) 


3.3 Consistency with Yeh’s theory for the small-period limit 


For the small-period limit the validity of these rules can be analytically verified by treating 
the periodic structure as alternation of two homogeneous layers, where we use the boundary 
conditions for the continuity of the tangential electric and magnetic fields on all interfaces. 
Now the function £ is assumed constant (er or £r) within each layer of the thickness dı or d2. 
According to the geometry in Fig. 1, the field in the jth layer has the dependence 


E x,y) = e—ikopox [q e—ikon;(y—y;) aie betkoui(¥—¥i)), (37) 


where qj = (£j — pa)/2, with €; = £q for odd j and £j = £q for even j. It is coupled with the 
field in the next layer by the matrix equation 


Aji. _ P; 0 a; B a; Bj 
salola 0) o- (RS) » 


J 
where a; = 5(1 7 qj/9j+41)» By = te — qj/qj+1) for the TE polarization, or a; = 5(1 + 


Ejlj+1/Ej+14j) Bj = 5(1— €;9j+1/€)+19;) for the TM polarization, and P; = e~ t04;4; for both 


1/2 


polarizations. Obviously, qj = (€j — p3)'/, assuming the e~‘*oPo* factor of all fields. 


Applying the small-period approximation Ce a ~ 1 + ikgqjd; to the problem of propagation 
through the whole period, 


a3] fa _ [RO P 0 
slefa) e=e[eala[e a]. g 


2 l 


yields the eigenvalues of the Q operator 


Desti ikoA/ £0 -p?, (TE) (40) 


Os = 1 ikyAy/(1— pBno)e0. (TM) (41) 
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with 
£o = wep +(1—w)ez, (42) 
1o = wyr + (1—w)nn, (43) 


Assuming a Floquet mode with the eigenvalues of Q being Q4 = eto04 ~ 1+ ikgqoA, we 
see that in the small-period limit the periodic structure behaves as a homogeneous anisotropic 
medium with the y-component of the normalized wave vector qo = (£0 — be) 2 for the TE 


polarization and qo = [(1 — p3yo)eo]!/? for the TM polarization. These formulae are identical 
with Equations 35 and 36 if we retain only the Oth element of all matrices. This is a very 
interesting disclosure that the results obtained by Yeh and coauthors already in 1970s are 
consistent with the extensive, more general research carried out in 1990. 


3.4 Numerical example 


Example of the comparison of applying the correct Fourier factorization rules with applying 
the opposite ones is shown in Fig. 2 for both polarizations. The normalized eigenfrequency 
wA/2rc of the first band is displayed according to N; the structure is made as two alternating 
layers of the equal thicknesses 500 nm (A = 1000 nm) and permittivities e; = 3 and ey = 1. 
The wave vector is chosen k = (0.57t/A){1,1, 0]. 


0.2472 
0.247 % —- correct factorization 9000000 OOF 
x -x- opposite factorization o eaii J 
> 0.2468} A 
c X 
D ae 
> Xx 
© 0.2466 x 
D la 
N i —-correct factorization 
g 0.2464 x -x- opposite factorization 
g x 
= i 
0.2462 i 
; 
0.246 ' 
J 
0 10 20 E 30 10 20 N 30 40 50 
(a) TE polarization (b) TM polarization 


Fig. 2. Convergence properties of the correct factorization compared with the opposite one 


4. Two-dimensional photonic crystals with circular elements 
4.1 Geometry of the problem 


Fig. 3 displays the geometrical arrangement of a 2D photonic crystal made as bi-periodic 
alternation of rods or holes with a cylindrical cross-section. Instead of a single vector of 
periodicity we now have the plane of periodicity (determined by two vectors of periodicities 
defining a unit cell), which here coincides with the xy plane. For simplicity we choose the 
incidence direction in the plane of periodicity, so that we can again distinguish between two 
independent polarizations. The plane of incidence is now determined along the propagation 
wave vector k = ko|po,qo,0| and perpendicular to the plane of periodicity. The TE 
polarization has now H along the z axis, which is now the more difficult case (with nonzero 
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plane of incidence 


swks 


"Mt nti a i le 


Fig. 3. Geometry of a 2D photonic crystal made as bi-periodic alternation of rods or holes 


Hz Ex, Ey), while the TM polarization has E along the z axis (nonzero Ez, Hy, Ay), which is 
now the simple case. 


0000o 0000o 
0000o "YYYY, 
00000 © 000 

Alleleceeee ~l eee 

, ilele 1@6@ 000 
—— 

P 
A penodie pal A periodic cell 

Reciprocal Reciprocal  % 

space: space: 


(a) Square periodicity (b) Hexagonal periodicity 


Fig. 4. Two examples of 2D periodic arrangements; below are the first Brillouin zones in the 
reciprocal space 


As an example, in Fig. 3 we have chosen the xz plane as the plane of incidence, but we 
could choose any plane parallel with the z axis provided that we want to treat the TE and 
TM polarizations independently. In general, the plane of incidence is determined by the k 
vector in the reciprocal space as displayed in Fig. 4, whose particular symmetry points are 
denoted T, X, M, or K according to the corresponding periodicity. 


In this section we study a 2D photonic crystal composed of infinite cylinders with a circular 
cross-section with either square [Fig. 4(a)] or hexagonal [Fig. 4(b)] periodicity. For the square 
symmetry the unit cell has the dimensions A, = Ay = A, and for the hexagonal symmetry 
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it can be chosen as an A-by-AV/3 rectangle. The corresponding first Brillouin zones of the 
reciprocal space are depicted in the bottom part of Fig. 4. 


Maxwel’s equation are again 
dyEz = —iwyugoHz, dx Ez = iwugoHy, Ox Hy — Oy Hx = iwegeEz, (44) 


However, now we cannot put there the labels (£) or (T) for the factorization type as easily 
as above, because the discontinuities of the permittivity function are now mixed among the 
components of E. 


Assuming a hypothetical anisotropy of the relative permittivity function, we define a scaled 
electrical displacement D, 


Dy, Ex Exx Exy | | Ex | 
~ = = P TE 46 
[ø e [E ie Eyy | | Ey i ea 
Ď- == €z7E- = E€) Ez, (TM) (47) 


where €zz is obviously the only component of the permittivity tensor for which we can use the 
Laurent rule. 


Defining also a 2-by-2 matrix of electrical impermittivity 7 = e7! helps in the formulation of 
the wave equations 


2 
— (32 +32) Ez = SewEz (TM) (49) 
where 77;x are the components of the electrical impermittivity. For the simplicity of the TM 
polarization case we below focus our attention only to the TE polarization. 
4.2 Methods of Fourier factorization 


In this section we compare several models corresponding to different factorization 
approaches. 


4.2.1 Elementary (Cartesian) method (Model A) 


First, Model A assumes the solution in the basis of the ĉ and g polarizations uniform within 
the periodic cell, where in accordance with Ho et al. (1990) we choose the Laurent rules 

[Dx] = [Ex] = [e] [Ex], (50) 
[Dy] = [eE] = [le] [Ey]. (51) 
The components of the electric impermittivity in Equation 48 then becomes [e] ~t} for the cases 
of 7xx, Nyy, and zero for the cases of xy, Nyx, OF 


[el~* [o] | 

n]a = ee (52) 
ma = [f fa 

For illustration we show the distribution of the first basis polarization vector (identical with 


the constant vector ĉ) in Fig. 5(a), where the black circle denotes the element boundary (the 
permittivity discontinuity). 


www.intechopen.com 


Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 329 


4.2.2 Normal vector method (Model B) 
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Fig. 5. Distribution of the basis polarization vector u for the factorization models 


According to the factorization theorems, neither the Laurent rule nor the inverse rule is correct 
for both products in Equations 50 and 51, because both pairs of functions have concurrent 
discontinuities and both products Ď, and Dy are discontinuous as well. On the other hand, 
by an appropriate change of the polarization bases at all points (using a space-dependent 


Jones matrix transform F), 
Ex Eu 
=F ; 53 
| | Ev | i ) 


we can treat independently the normal (u) and tangential (v) components of the fields by the 
correct rules, 


[Du] = [17e] 7+ [En], (54) 
[D>] = le] [Ev]. (55) 
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The field components E,,, D, are normal to the discontinuities of the relative permittivity 
function, while E,, Dy are tangential. The factorization rules used in Equations 54 and 55 are 
justified simply because E, and D,, are continuous. 


A suitable distribution of the matrix F within the periodic cell can obviously be the rotation 


F= Epo (56) 


sing coso 


where the polar angle (x, y) is in the first cell distributed according to the polar coordinates 
re’? = x + iy, and then periodically repeated over the entire 2D space. This enables defining 
the matrices |c], |s] from the corresponding 2D-periodic functions c = cos ¢, s = sin Q. 


Let u and v be the two columns of the matrix F, both being mutually orthogonal basis vectors 
of linear polarization. From the above definitions we see that u is a polarization vector normal 
to the structure discontinuities, whereas v is tangential. In Fig. 5(b) we show the distribution 
of u within the periodic cell. The basis polarization vectors are constant along the lines of the 
constant azimuth (¢ = const) and rotate as ¢ increases. It is obvious that the matrix function 
F(x,y) has no discontinuities concurrent with the electric field, so that we can use both Laurent 
and inverse rules for the transformation of polarization, e.g., 


fel] =m [et i 


_ [fel ts] 
m= [iy a | a 


Combining Equations 54, 55, and 57 yields 


al p [0] | “a | 

= |F] | £ F a | a 59 
[te = (iaf fe | oy Pa 
from where we derive the electric impermittivity in the reciprocal space (corresponding to 
Model B) 


tole = 11 [G L, | oe 


= vie + [el "Is". [elles] - iia 
Lelles] — lel ies], Lells] + fel~ eT] ’ 


whose components are immediately applicable to Equation 48. 


(60) 


4.2.3 Method with elliptical polarization bases (Model C) 


The above approach (Model B) only deals with linear polarizations and thus suffers from 
the fact that the matrix function F(x,y) has a singularity at the central point of the periodic 
cell and other discontinuities along the cell boundaries. This slows down the convergence of 
the numerical implementation, as will be evidenced below. On the other hand, we can make F 
continuous by using complex functions ¢ and Ç or, in other words, by defining u, v as complex 
vectors corresponding to generally elliptic polarizations, 


e a 
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(still orthogonal). By means of rotation @ and ellipticity e we define the first basis vector 


ig | cos@ —siné COS € 
aT eae cos 6 | bal á (62) 
where 
A(r, ) = ~, (63) 
7T Tr 
\ % (1+cos ¥) (r< R) 
elr,o) c {1+ cos TERR (r >R). (64) 


Here R denotes the radius of the circular element and 


A/2 
max(| cos $|, | sin |) 


D(p) = (65) 
is the distance from the cell’s center to its edge. In Equation 62 the Jones vector on the right 
represents a polarization ellipse (with ellipticity €) oriented along the x coordinate, the matrix 
in the middle rotates this polarization by the azimuth 6, and the factor e® preserves the 
continuity of the phase at the center and along the boundaries of the cell. This continuity 
can be easily checked by evaluating the limits 


1/1 
limu= li PEFEA y 
r0 De VA Bg a 


which is the vector of left circular polarization (independent of ¢). 


The distribution of the basis polarization vector u within the periodic cell is shown in Fig. 5(c). 
Here the azimuth of the polarization ellipse is constant along the lines coming from the cell’s 
center, which is similar to Model B. However, the ellipticity is now zero (corresponding to 
linear polarization) only on the boundaries of the circular element, has the maximum value 
(7/4 for circular polarization) at the cell’s center and along its boundaries, and continuously 
varies (with a smooth sine dependence) in the intermediate ranges. Thus we obtain a smooth 
and completely continuous matrix function F(x, y), which is analogously used to calculate the 
impermittivity in the reciprocal space 


inpe = EENEI + e1- ee"D, CADIE - fe Ee) a 
IMIE] [el “1661. ENII + bel ieg] 

In the case of the hexagonal periodicity we define u and the other periodic quantities inside 

one hexagon (half the area of the rectangular unit cell) where we can use formally the same 

equations as above, except for 


A/2 


_ my)’ 
max, [cos ($ — 3) ] 


D(¢) = (68) 


which is now the distance from the hexagon’s center to its edge. Here A is the hexagon’s 
shortest width (equal to the width A, of the rectangular cell). 
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4.2.4 Modified method for densely arranged elements (Model C’) 


To analyze a more complicated situation, we consider a photonic crystal with square 
periodicity where circular elements are densely arranged near each other, i.e., where the 
radius R is almost the half width A/2 of the periodic cell. Then the convergence properties of 
F becomes worse, which affects all the derived quantities. For this reason we again redefine 
the polarization distribution. For the modified Model C’ we define u to be still same inside the 
circle (r < R), but different outside. Assuming the rotation and ellipticity along the boundary 
of the square cell 


(p) = 0 (D(¢),) = Fround (~/F), (69) 
Elp) = € (D(¢), p) = g (1 — cos4g) (70) 


(where “round” denotes rounding towards the nearest integer), we define the rotation and 
ellipticity outside the circle (r > R) as 


(1,0) = 4 {O6(h) +p + [6 (p) — p] cos TRER, (71) 
e(r,p) = slg) {1 + cos ee ra \ ; (72) 


Assuming otherwise the same Equations 59, 61, 62, and 65, we obtain for |n] c formally the 
same matrix as in Equation 67, except that the functions ¢ and ¢ are now derived from different 
azimuth and ellipticity distributions of u. Note that u is again continuous along the cell’s 
boundaries; to evaluate its precise limits [when x => +D(0), y = const or y > +D(¢/2), 
x = const] would now be more complicated. The distribution of the basis polarization vector 
u within the periodic cell is depicted in Fig. 5(d) together with dimensions. 


4.3 Numerical examples 


We examine the numerical performances of all factorization models presented in Section 4.2 on 
three samples of 2D photonic crystals, for which we calculate the eigenfrequencies wx (where 
the band number x = 1 stands for the lowest eigenfrequency, x = 2 for the second lowest, etc.) 
and the corresponding eigenvectors |Hz]x of Equation 48. All convergences will be presented 
according to the maximum Fourier harmonics retained inside the periodic medium, which 
will be kept same for the x and y directions (M = N). 


First, Sample S1 is a square array of cylindrical rods of the circular cross-section with the 
diameter 2R = 500 nm, square period A = 1000 nm, relative permittivity of the rods £1 = 9, 
and relative permittivity of the surrounding medium corresponding to vacuum (e2 = 1). Its 
dispersion relation is displayed in Fig. 6(a). Similarly, for Sample S2 we assume exactly the 
same parameters except the diameter of the rods, now being 2R = 900 nm. This corresponds 
to densely arranged elements (the distance between two adjacent rods is only 100 nm). Finally, 
for Sample H we consider a hexagonal array of cylindrical holes of the circular cross-section 
with the diameter 2R = 600 nm, hexagonal periodicity A = 1000 nm (corresponding to the 
rectangular cell of the dimensions Ay = 1 um, Ay = V3 um), relative permittivity of the 
holes corresponding to vacuum (£1 = 1), and relative permittivity of the substrate medium 
(surrounding holes) £2 = 12. The dispersion relation of Sample H is displayed in Fig. 6(b). 


For our analysis we choose the eigenmodes T2 and T3 of Sample S1, the eigenmode X3 of 
Sample S2, and the eigenmode M1 of Sample H, where the letter denotes a point of symmetry 
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Fig. 6. Photonic band structures of two samples 


Fig. 7. Amplitude distributions of the scaled magnetic field |Hz| within the unit cell of four 
chosen eigenmodes 


and the number corresponds to the band. The amplitude distributions of the modes are 
displayed in Fig. 7, and the corresponding convergence properties of the eigenfrequencies 
calculated by the above described models are shown in Fig. 8 (Model C’ is only compared for 
Sample S2). 


The result of a careful comparison of the numerical efficiencies of all the factorization models 
can be summarized as follows. Models B, C, and C’ always converge considerably faster 
than Model A. Model C converges faster (having usually one order higher precision) than 
Model B with two exceptions. The first exception is the case such as in Fig. 7 (S1-T2) and 
Fig. 8(a), where the discontinuities of the polarization transformation matrix F coincide with 
a nearly zero amplitude of the field, so that the discontinuities do not manifest themselves. 
The second exception is the case such as Fig. 7 (S2-X3) and Fig. 8(d), where the elements 
are densely arranged which causes rapid variations of the ellipticity between two adjacent 
elements (which are very close to each other); this requires more Fourier components than the 
weak discontinuity of the linear polarization u in Model B. The problem is solved in Model C’, 
which obviously converges fastest among all the four models applied to Sample $2. 


www.intechopen.com 


334 Photonic Crystals — Introduction, Applications and Theory 


0.52582 I sess, 
xx | 
0.52581 0.7188 f 
i x- model A 

3 o5258} | ~~ model B a l 
Q —— model C g H 
S S 0.7184 Í 
g 0.52579 3 i 
5 5 
8 0.52578 g iN Y 
© æ 0.718 x- model A 
E 0.52578 E ~-a- model B 
= = —— model C 


0.52577 
0.7176 


0.52576 


(b) Sample S1, point T3 


0.1842 p 0.4422 


0.1841 L 


> > 
(2) O 
(= H Cc 
3 S 
2 0.184 | > 
£ l © 
Ba > 0.4421 
Oo : oO 
DN 3 N 
= 0.1839 = 
£ £ 
g g 
iy x- model A --&-- model B 
0.1838 h = model B —— model C 
—— model C —— model C' 
0.1837 0.442 Hi 
5 10 15 20 25 30 
N 
(c) Sample H, point M1 (d) Sample S2, point X3 


Fig. 8. Convergences of the Fourier factorization methods for three samples 


5. Two-dimensional photonic crystals with non-circular elements 


In this section we briefly show numerical efficiencies of three factorization models derived 
above (Models A, B, and C) applied to 2D photonic crystals made as long elements of other 
shapes, namely rods with the square cross-section and tubes with the ring and split-ring 
cross-sections, arranged with the square periodicity. For all the three samples we choose the 
permittivity of the elements € = 9 and permittivity of vacuum ¢ = 1 for the surrounding 
medium. 


5.1 Periodic rods with the square cross-section 


For the photonic crystals made of square-sectioned rods, the period is chosen A = 1000 nm, 
and the width of the square d = 600 nm. The distribution of the basis polarization vector u, 
analogously to Section 4, is displayed in Fig. 9 for all the three compared models. 


For Model A, of course, u = ĉ as visible in Fig. 9(a). For Model B, we divide the unit cell into 
four parts by its diagonals (the lines x = y and x = —y). The vector u(¢), depending only 
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Fig. 9. Polarization distribution of the factorization methods for periodic squares 


on the azimuthal angle, is set ĉ for ¢ € (— F, 4) U (32 ; Sz ), and g for the remaining angles. 
This means that u is always linear, perpendicular to the permittivity discontinuities, and its 
discontinuities (along the lines x = y and x = —y) have no concurrent discontinuities with 
the permittivity function except those at the four corners of the square. Hence, Model B here 
fulfills nearly the same conditions for the application of the factorization rules as demanded 
in Section 4 for circular elements. 


For Model C we divide the unit cell into four areas in the same manner. This time, however, 
the basis vector u(r, p) depends on both polar coordinates and the corresponding polarization 
is in general elliptic. In analogy with Section 4 we want u to be perpendicular to the 
permittivity discontinuities and to remove its discontinuities as much as possible. The most 
simple way how to do this, although the discontinuities will not be completely removed, is 
to set the azimuth 6(#) of the polarization to zero for p € (—4, 4) U (22, %2), and F for 
the remaining angles, and to use Equation 64 for the ellipticity distribution, where D(p) now 
corresponds to the square element. 
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Fig. 10. Convergence of the factorization methods for rods with the square cross-section 
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The numerical efficiencies of all the three models are displayed in Fig. 10 for the 3rd band of 
the symmetry point T. As clearly visible, Models B and C converge considerably faster than 
Model A, with Model C being a little better. Although the improvement from Model B towards 
Model C is not so distinct, it could be improved by a more careful choice of a completely 
continuous distribution of the polarization basis. 


5.2 Periodic hollow cylinders 


For the photonic crystals made of periodic hollow cylinders (tubes with the symmetric ring 
cross-section) we choose the period A = 1000 nm, the inner diameter (the diameter of the 
inner circular hole) R4 = 400 nm, and the outer diameter R} = 680 nm. The distribution of 
the basis polarization vector u for all the three compared models is displayed in Fig. 11. 


(a) Model A 


Fig. 11. Polarization distribution of the factorization methods for hollow cylinders 


For Model A again u = ĉ. For both Models B and C we use the same polarization distributions 
as in Section 4 for the inner area (r < R1) and for the outer area (r > R>). For the annulus 
area (Ry < r < R2) we simply choose the polarization distribution same as that in Model B. 
These distributions are obviously the most straightforward analogies of the distributions used 
in Section 4 for circular elements. 


The numerical efficiencies of all the three models are displayed in Fig. 12 for the 3rd band 
of the symmetry point T. As clearly visible, Models B and C converge considerably faster 
than Model A, but now Model C exhibits no improvement against Model B. This is because 
the discontinuities of w at the center and along the boundaries of the periodic cell quite well 
coincide with the zero amplitude of the mode, as visible in the inset of Fig. 12, so that the 
discontinuities do not manifest themselves in the calculations. 


5.3 Periodic split hollow cylinders 


For the photonic crystals made of split hollow cylinders (tubes with an asymmetric, split ring 
cross-section) we choose the period A = 1000 nm, the inner diameter (the diameter of the 
inner semi-circular hole) R = 600 nm, the outer diameter Rə = 720 nm, and the relative 
azimuthal length of the ring wọ = 0.9 (where wp = 1 means the complete, symmetric ring). 
The distributions of the basis polarization vector u for all the three compared models are 
displayed in Fig. 13. 
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Fig. 13. Polarization distribution of the factorization methods for split hollow cylinders 


Since the area of splitting is quite small compared to the full area of the unit cell, we have 
chosen the polarization distributions of all the three models exactly same as for the symmetric 
rings in Section 5.2. The condition for normal u and tangential v is not satisfied on the surface 
proportional to the length 2 x 120 = 240 nm, but is satisfied on the surface proportional 
to 0.9 x 27æ(R1 + R2) % 7464 nm, which is 30 times higher area, justifying this negligence. 
Of course, for modes with most of electromagnetic field resonating in the critical area this 
approximation would be insufficient. 


The numerical efficiencies of all the three models are displayed in Fig. 14, again for the 3rd 
band of the symmetry point T. We can describe these performances by exactly the same 
conclusion as for the symmetric rings in Section 5.2. Models B and C converge similarly and 
both considerably faster than Model A. The discontinuities of u at the center and along the 
boundaries of the periodic cell again coincide with the zero amplitude of the mode, though 
with some deviations visible in the inset of Fig. 14. 
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Fig. 14. Convergence of the factorization methods for split hollow cylinders 


6. Anisotropic photonic crystals 


In this section we will briefly demonstrate the application of the factorization rules to 2D 
photonic crystals made of anisotropic materials, again with the plane of incidence parallel 
to the z axis (with geometry of Fig. 3). Unlike the isotropic crystals, now the TE and TM 
polarizations are not separable. Instead of the scalar permittivity we define and expand the 
components of the relative permittivity tensor function 


00 
ae = ) ener ™, (73) 


mn=—oo 


where €jkmn are the Fourier coefficients. The wave equation for a generally anisotropic 
medium, now described by the permittivity tensor 


Exx Exy Exz 
E = | Eyx Eyy Eyz |, (74) 
Ezy Ezy Ezz 
is the operator equation 
>20; Ox0y 0 
CE = E Ĉ=e! | 0x0, -Ə 0 . (75) 


0 0 -32-a 


Similarly as above, Model A assumes all components of the permittivity tensor treated by the 


Laurent rule, i.e., . 
[Dj] = $ Jex] [Ez]: (76) 
k 
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To apply the factorization correctly, we must again separate the normal and tangential 
components of all fields for which we can use the correct rules. Let us define a 
space-dependent matrix transform F(x, y) so that 


Ex Eu 
Ez Ez 


where E, and E, are the normal and tangential components of the vector E to all 
discontinuities of the permittivity. 


Analogously to Section 4, for the case of circular elements we can choose the polarization basis 
distribution corresponding to Models B and C as 


cosh —sing 0 
F= |sinọġ cosg 0|, (ModelB) (78) 
0 0 1 
= 0 
F= | ¢* 0}, (Model C) (79) 
0 0 1 


where the matrix elements have the same meaning as in Section 4. In the new coordinates we 
can write 


Dd, Euu Euv Euz Ey 
Dy = | €ou Evo Evz Ev |. (80) 
Dz Ez 


Ezu Ezv Ezz 
Now let us separate two sets of quantities, those which are continuous (D,,, Ev, Ez) and those 
which are not continuous (Dy, Dz, E,,) to the discontinuities of the permittivity. Expressing 
the second set according to the first one yields 


1 — Ew — Euz 
x Euu Euu Euu 
Eu Du 
D È EvuE EvuE 
Dy = G E; ; G + Evy a uv. Ez z uz E (81) 
dD E uu uu uu 
zZ PA 
Ezu Elo EzuEuv Ez EzuEuz 


Euu Euu 


for which we can simply use the Laurent rule, 


[Eu] [Du] 
[Do] | = [G] | [Eo] |- (82) 
[Dz] [Ez 


From this we express [Dj] according to [E;], 


[Du] [Eu] 
[Do] | = le{uoz3] gc | [Eo] | > (83) 
[Dz] [Ez] 
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with 


a 1] fe 
Euu á Euu Euu ||” 
—1 —1 
Evu na € EvuEuv Evu 1 Euv 
Euu Euu d ue Euu Euu Euu Euu 
T 


$ 
y 


+l 
= E EzuEuv + Ezu 1 l Euv 
£ y Euu Euu Euu Euu 
A || Euz 
Euu Euu 
-1 
EvuEuz Evu 1 Euz : 84 
[eve Euu | + [ Euu | [2 | I Euu | ( ) 
EzE E 1] ‘Te 
ZUSUZ Zu Suz 
lezz Euu | + | = | È | | | 
Finally we obtain the matrix of permittivity in the Cartesian coordinates via the formula 


le geyz}]p c = Mlle uoz} lg, cE] (85) 


where the index B or C corresponds to the chosen model. 


7. Distribution of polarization basis for 3D structures 


In this section we will suggest how to create the polarization basis distribution for a 3D 
photonic crystal. Unlike the previous cases, we now have a fully vectorial wave Equation 1. 
The corresponding material equation D = e(x, y, z) E must be changed to separate the normal 
and tangential components of the fields to the e discontinuities, which are now surfaces in 
the 3D space. For simplicity we assume a 3D photonic crystal made as spheres (with the 
radius R) arranged in the space with the cubic periodicity. To make a 3D analogy with Model C 
described in the previous sections, we must find a matrix transform F(x, y, z) whose columns, 
denoted u, v, and w, are complex vectorial functions of space, mutually orthonormal and 
continuous at all points, where u is the normal vector at each point of the sphere’s surface and 
v and w are tangential. 


As the first step we choose the distribution of these vectors on the sphere. Defining the 
spherical coordinates 


x = rsin cos ®@, (86) 
=rsinvsing, (87) 
z=reosd (88) 


A 


and the corresponding unit vectors 7, 3, and @ (pointing along the increase of the 
corresponding coordinate) helps us to define the vectors 


ur(8,p) = u(R, 0,9) = Ft f, (89) 


; 1 n A 
— — Otep — S o : 

vR(0,¢) = v(R,0,¢) =e Laor e (90) 

1 


V, — R, V, = i(O+p) 
we (8,9) = w(R,8,p) = 0+#) c 


(3 cos 0 — id). (91) 
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Obviously, u is normal on the whole sphere’s surface. The vector v corresponds to the left 
circular polarization at the sphere’s poles (z = +R), to the vertical linear polarization on the 
sphere’s equator (z = 0), and is elliptical and continuous in the intermediate ranges. The 
vector w is simply chosen orthogonal to both u and v. Both v and w are tangential to the 
sphere’s surface. 


If we extend the radial dependence of the vectors defined by Equations 89-91 to the entire 
cubic cell (formally replacing R by r), then we obtain some analogy with Model B. The 
obtained matrix F has no concurrent discontinuities with g, but it is discontinuous at the center 
and on the boundaries of the cubic cell. To remove the discontinuity at the cell’s center, we 
will proceed differently. 


As the second step, we define the basis vectors at the center of the sphere, 


uo = u(0,8,9) = -5 (8 — i8), (92) 
vo = v(0,8,ġ) = 2, (93) 
wy = w(0,8,9) = (ê + ið), (94) 


which are again mutually orthogonal. Then we extend u onto the whole cubic cell, 


(uo + ur) + 4 (uo — ug) cos (r< R) 


= 95 
(uo + ur) + £ (uo — ug) cos EPR) (R << D). 9) 


1 
dhext (T, v, op) = i A 
2 


where D(0,¢) is the distance from the cell’s center to its boundary along the ray determined 
by the spherical angles. The desired unit vector of the polarization basis is then 


u = AwUext, (96) 


where Au = [4 (1 + cos? 4)]~1/? for r < Rand Ay = [3(1+ cos? ae forr > R 
is a scalar function ensuring that u becomes a unit vector everywhere. We could extend v and 
w in a similar way, 


Í (vo + vr) + (vo — vg) cos & (r< R) 
dgis 7\VO 200 R 7 a 97 
ie | 3 (vo + vr) + $ (vo — vR) cos AR (R<r<D). a 

5 (wo + wr) + 4 (wo — wR) cos ¥ (r< R) 
rla) | I\%0 2 A T 98 
Wext( $) l 5 (wo + wr) + 5(wo — wR) co ae ee (RTS DP). i 


but it is not clear whether vext OF Wext are both perpendicular to u and mutually. To ensure it 
we define 
v = Ay(1 m Pu) Vext, (99) 


where Pu = uu is the projector into the space of vectors proportional to u, so that 1 — Pu is 
the projector to its orthogonal complement. Similarly, 


w = Aw(1 =Pu= Py) Wext, (100) 


where 1 — Pu — P, is the projector into the space of vectors perpendicular to both u and v. 
Here A, and Aw are analogously chosen scalar functions ensuring the unit size of the 
corresponding vectors. 
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Finally, the matrix of permittivity in the reciprocal space with correct Fourier factorization, 
directly applicable to the wave Equation 1, becomes 


(1° (0) fol] 
[elc = [F] | fo} [e] lo] | IF 1. (101) 
fo} fo) fel 


because the first element on the diagonal corresponds to the normal (u) components of the 
fields and the other two correspond to the tangential components (v, w) of the fields. 


8. Conclusion 


We have derived the methodology how to apply the Fourier factorization rules of Li (1996) to 
various photonic crystals. For the case of 1D crystals there is clear consistency of the correct 
rules with the classical theory of Yeh et al. (1977). For 2D crystals the convergence properties 
strongly depend on the chosen distribution of the polarization basis; we have shown that it 
is desirable to choose a distribution as smooth as possible. The method is also usable for 
periodic elements of any shape, where complicated shapes require complicated distributions 
of polarization bases. We can also use it to simulate 2D periodic elements made of anisotropic 
materials, as well as 3D periodic crystals. Moreover, the method can also be used to photonic 
devices such as photonic crystal waveguides by applying the demonstrated methodology 
to the device supercell. It is particularly advantageous for devices where high accuracy is 
required, e.g., for analyzing defect modes near photonic band edges [Dossou et al. (2009); 
Mahmoodian et al. (2009)], and for large devices for which the available computer memory 
enables calculations with only a few Fourier components (photonic crystals fibers with large 
cladding, or asymmetric 3D crystals). 
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